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The question of whether electron-electron interactions can drive a metal to insulator transition in 
graphene under realistic experimental conditions is addressed. Using three representative methods to 
calculate the effective long-range Coulomb interaction between 7r-electrons in graphene and solving 
for the ground state using quantum Monte Carlo methods, we argue that without strain, graphene 
remains metallic and changing the substrate from Si 02 to suspended samples hardly makes any 
difference. In contrast, applying a rather large - but experimentally realistic - uniform and isotropic 
strain of about 15% seems to be a promising route to making graphene an antiferromagnetic Mott 
insulator. 
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Over the past decade graphene has established itself 
as a remarkable new material with superlative properties 
HIE]. However, the early hopes to utilize it as a next 
generation transistor have been dashed mostly because 
graphene remains metallic - these prototypical Dirac 
fermions are immune to many of the conventional routes 
for driving two-dimensional electron gases into an insu¬ 
lating state, including, for example, Anderson localiza¬ 
tion and percolation transitions (see e.g. Ref. 0). Other 
mechanisms for opening band-gaps including hydrogena¬ 
tion |4], application of uniaxial strain [5] and forming 
nanoribbons [6] severely degrade graphene’s mobility. 
Very recently, moire heterostuctures using graphene and 
hexagonal boron nitride have shown evidence of an in¬ 
sulating phase M, although interpreting these results 
remains somewhat controversial pus. 

In this Letter, we explore a different avenue to 
make graphene insulating, namely, utilizing the electron- 
electron interactions. Despite much study on the ef¬ 
fects of interactions in graphene m it is surprising how 
much still remains to be understood. While it is clear 
that without any electron-electron interactions, graphene 
should be a semi-metal (SM), and that for very strong 
interactions it should be an insulating anti-ferromagnet 
(AFM), it remains unclear what one should expect for the 
real graphene material. For example, there are distinct 
claims in the literature that suspended graphene should 
be insulating, strongly metallic and weakly metallic [H- 
16 . This discussion could have practical relevance as it 
could be the basis for a low power Mott-transistor m- 

In this work we explore different ways of controlling the 
effective strength of electron-electron interactions in re¬ 
alistic graphene devices, and propose how one can move 
around its phase diagram. In particular (and in con¬ 
trast to what is widely assumed to be true mm) , we 


demonstrate that it is the non-universal, material-specific 
and short-range part of the electron-electron interactions 
that plays the dominant role in determining graphene’s 
ground state. More interestingly, we conclude that ap¬ 
plication of isotropic strain is considerably more efficient 
in approaching the SM-AFM phase transition than sub¬ 
strate manipulation, providing a new route for driving 
the system into the elusive Mott insulating phase that 
has yet to be observed experimentally. 

The Hubbard model has served as a versatile paradigm 
to study interacting electrons on a lattice. It is defined as 
an effective model for electrons in partially filled narrow 
energy bands of a crystal’s spectrum. While the canoni¬ 
cal Hubbard model includes only on-site interactions, the 
effects of longer range interactions are incorporated by a 
straightforward generalization of the two-body interac¬ 
tion term, described by the Hamiltonian 

H = -t Y + hx -) + F 

<U),cr i 

+ 2 y! niaVijnja' 5 (1) 

where cj a (q ct ) creates (annihilates) an electron of spin 
a =tl at position while hi a = cj a Q a gives the density 
of electrons with spin a at position r*. The nearest neigh¬ 
bor hopping integral is identified by £, while Vij stands 
for the interaction between electrons at sites i and j. We 
note that a realistic description of graphene requires the 
parameters Vij to be fixed in accordance with the spatial 
profile of the (partially screened) Coulomb interaction 
V(r) that results from the screening of the bare Coulomb 
interaction by electrons in energy bands other than the 
7r-bands. 

It is well established that the canonical Hubbard model 
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FIG. 1. Schematics of biaxially strained graphene (left panel). 
Representation of graphene’s low-energy spectrum (right pan¬ 
els) of unstrained non-interacting and interacting graphene at 
half-filling, and biaxially strained non-interacting and inter¬ 
acting graphene at half-filling. 


[Eq. (§ with only on-site interactions, i.e., Vu = U (> 0) 
and Vij = 0 for all i % j] on the honeycomb lattice at 
half-filling has a critical Hubbard on-site interaction pa¬ 
rameter for the SM-AFM transition U c = (3.80 =b 0.01)£ 
jT8H20l . Following various works based on ab initio meth¬ 
ods (see e.g. Ref. [21] [22]). it is generally agreed that 
t = (2.7 ± 0.2) eV. Estimates of the on-site interac¬ 
tion parameter U for realistic experimental realizations 
of graphene vary widely in the literature [231126] . with 
values ranging from U ~ 1 eV to 10 eV (where the lower 
estimates would suggest that graphene is metallic, while 
the higher estimates hint at it being insulating). 

However, ignoring longer range interactions in 
graphene is problematic since the long-range tails of the 
Coulomb potential between Dirac fermions cannot be ef¬ 
ficiently screened [27] . To address these Coulomb tails, 
it was recently conjectured m that the effects of non¬ 
local interactions can be mapped into the Hubbard model 
with an on-site interaction U given by U ~ U — V, where 
U = Vu corresponds to the on-site interaction of the 
long range Hubbard model, while V = Va +5 stands for 
the value of the Coulomb potential between electrons at 
two neighboring atoms on the honeycomb lattice. This 
effective U would thus be the crucial factor determining 
graphene’s phase. As we discuss below, our numerical 
calculation with the full long-range potential shows that 
this approximation is qualitatively correct, but quantita¬ 
tively inaccurate. 

Here we study the possibility to drive graphene across 
the SM-AFM phase transition by substrate manipula¬ 
tion or application of biaxial (i. e. uniform and isotropic) 
strain - see Fig. [I]). First we must fix the long range 
Hubbard model’s parameters Vij. These are the crucial 
ingredients determining the ground state properties of 
the system, yet their real values are unknown. We use 
three representative methods to capture the full spatial 


profile of the partially screened Coulomb interaction for 
p z electrons in realistic graphene, and choose Vij accord¬ 
ingly. These methods will be discussed in detail below, 
but now we just introduce their names: Thomas-Fermi 
(TF), constrained random phase approximation (cRPA) 
and the quantum chemistry - Pariser-Parr-Pople (QC- 
PPP) method. We then investigate the effect of biax¬ 
ial strain and substrate manipulation on the partially 
screened Coulomb potential V(r). We find (see Fig. [2| 
that biaxial strain strongly modifies the V(r) close to 
r = 0 (not affecting the long-range interactions), while 
changing the substrate (which changes both the dielectric 
screening [28] and the amount of disorder 0 ) only weakly 
modifies the long-range tail of V(r). Finally, using quan¬ 
tum Monte Carlo techniques (finite temperature deter¬ 
minant quantum Monte-Carlo and zero-temperature pro¬ 
jective quantum Monte-Carlo), we simulate the ground 
state of the long-range half-filled Hubbard model (in the 
honeycomb lattice) with the Vij obtained from V(r), and 
argue that at least within the Thomas-Fermi approxima¬ 
tion, an experimentally feasible [29, 30. amount of strain 
would drive graphene into an interaction driven insulat¬ 
ing phase, which could be then measured in compress¬ 
ibility, transport or scanning probe experiments. 

We now detail the three methods that we use to es¬ 
timate the partially screened Coulomb interactions that 
p z electrons feel. (These were chosen since they are very 
representative of the different approaches that have so 
far been used in the literature). The cRPA method (see 
e.g. Ref. [31] ) was adapted to graphene by Wehling 
et al. [25]. In a systematic way, this method makes 
use of the electronic structure of graphene to compute, 
within the Random Phase Approximation, the polar¬ 
ization function P r associated with all the interaction 
events other than those involving two electrons from 
the 7r-bands. Then, the effective (partially) screened 
Coulomb interaction felt by the p z electrons is given by 
V(r) = [1 - V b are(r)P r \~ 1 V b are{r), where Vbare(r ) stands 
for the bare Coulomb potential. The accuracy of this 
method has long been debated in the literature (see e.g. 
Ref. [32j), and its results are often difficult to interpret 
physically. For graphene, the fact that the Dirac band 
spans a broad energy window further complicates the ap¬ 
plication of the cRPA formalism. Notwithstanding these 
difficulties, the cRPA remains the best numerical tech¬ 
nique at our disposal to determine the Vij for graphene. 
In this manuscript we use the cRPA results previously 
obtained in Ref. [25] . which compute C7, V and t for bi¬ 
axial strains up to 12%. In this range of strains all these 
parameters show a linear behavior. In order to obtain 
the cRPA values of U,V and t for up to 18% strain (see 
Fig. §l) we have assumed that this linear behavior re¬ 
mains unchanged, extracting 17, V and t from a linear fit 
to Ref. [25]’s numerical results. 

The QC-PPP method (see e.g. Verges et al. [33] ) 
works by using ab initio Hartree-Fock and post-Hartree- 
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FIG. 2. Effect of biaxial strain (left panel) and substrate (right panel) on the partially screened Coulomb interaction. We 
use three representative models: constrained Random Phase Approximation (cRPA), circles/full curves; quantum chemistry - 
Pariser-Parr-Pople (QC-PPP), squares/dashed curves; and Thomas-Fermi (TF), triangles/dot-dashed curves, (a) Suspended 
graphene both unstrained and subject to 18% biaxial strain, (b) Unstrained graphene both suspended and deposited on SiCU 
compared to the bare Coulomb potential. 


Fock techniques to solve for the ground state energy of 
molecules comprising a small number of benzene rings. 
These energies are then compared to an exact diagonal- 
ization of the long range Hubbard model where the Ohno 
interpolation formula, V(r) = U/yj 1 + (Ur/e 2 ) 2 , is as¬ 
sumed for the Coulomb interaction. The V(0) = U is a 
free parameter that is fixed by requiring the minimization 
of the root-mean square of the ground state energy of the 
ab initio calculations and that of the long range Hubbard 
model. The QC-PPP values of V and V used in this 
manuscript were extracted from Ref. [33], which calcu¬ 
lates V(r) for the phenalenyl ( 3 H-C 13 H 9 ) molecule. This 
method gives an an upper bound for the Hubbard U in 
graphene since larger molecules would have more screen¬ 
ing and reduced V(r). Both the validity of the Ohno in¬ 
terpolation and the extrapolation to larger system sizes 
give some reasons for caution. It has nonetheless proven 
extremely useful for small 7r-conjugated planar polycyclic 
aromatic hydrocarbons comprising tens of atoms such as 
anthracene and polyacenes [33l [34] . 

Finally, inspired by the work of Jung and MacDon¬ 
ald [26] we have constructed a Thomas-Fermi model 
to account for the screening of higher energy bands in 
graphene. Within the Thomas-Fermi screening approxi¬ 
mation the on-site interaction U is given by 

p 2r p -k 0 \r 1 -r 2 \ 

u =l-J d 3 r!d 3 r 2 |*(n)| 2 |ri _ ra| |0(r 2 )| 2 , (2) 

while the Coulomb interaction between two 7 r-bands’ 
electrons positioned at neighboring atoms (distance S) 


V is given by 

p 2r -k 0 |ri—r 2 | 

v =4 d ^ ra ^ + 15)12 i^)i 2 • 

(3) 


Here, 0(r) stands for the p z -orbital’s wave-function 
(which we approximate by that of atomic hydrogen). The 
free parameter fco in Eqs. © and © is fixed by requiring 
that the hopping integral 


t = J d 3 r </>*(r + S) 


fi 2 V 2 

2rn 



e -fco|r-Rd " 

|r-R;| . 


<K r ), 
(4) 


is equal to the literature accepted value of to = 2.7 eV 
pH] . In parallel with what we do for the other two meth¬ 
ods, we then interpolate between V^’s short-range val¬ 
ues U and V and the long-range tail of Vij (see below). 
The procedure used to compute Vij of biaxially strained 
graphene is similar to that discussed earlier [34] . 

The computationally demanding method employed 
prevents us from simulating large size systems. In par¬ 
ticular, one must include the effect of the surrounding 
electrons since their inter-band polarizability contributes 
at all length scales | 2 ] thus modifying the effective di¬ 
electric constant from 1/r to l/[r(l + 7rr s /2)], where 
r s = 2 e 2 / [(fta+ttb)/^] is the effective fine structure con¬ 
stant (where K a and ^ are the dielectric constants above 
and below the graphene flake). The presence of disorder 
in the substrate can also be accounted for by introducing 
a modified screening function (see e.g. Ref. 23)- The 
full profile of the partially screened Coulomb interaction 
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is obtained by interpolating between the short-range re¬ 
sults , at first neighbor distance, and the long-range tail 
(assumed to start at the fourth neighbor). 

As we can see in Fig. [2j the short-range part of the 
partially screened Coulomb interactions V(r) is strongly 
affected by biaxial strain (left panel), while its long-range 
tails are nearly insensitive to strain. Manipulating the 
substrate has a very weak effect on the long-range tails of 
the partially screened Coulomb interactions (right panel). 

With the electron-electron interaction profiles of Fig. 
[2] we have fixed the long range Hubbard model’s parame¬ 
ters Vij and explored its ground state using auxiliary field 
quantum Monte Carlo simulations (made possible by re¬ 
cent works [IjBJ 36, 37]) - a numerically exact method 
for investigating strongly correlated systems. In this 
manuscript, we use different implementations of the aux¬ 
iliary field quantum Monte Carlo technique: finite tem¬ 
perature determinant quantum Monte-Carlo (DQMC), 
whose correlation functions are given by 

( 6 ) = lTr[def> H ] = ^ V[^ T \e- s ^d[^ T ],(b) 

(we refer the reader to Ref. [38] for details); and zero- 
temperature projective quantum Monte-Carlo (PQMC) 
(for details see e.g. Ref. [39]). where the correlation func¬ 
tions are given by 


~ (gojgjgo) <fr r |e-eA/ 2 o e -eg/ 2 |^ r ) 

W <$o|*o> e^°c <tf r |e-eq* r > 

In both cases we use a Hubbard-Stratonovich transforma¬ 
tion to convert the interaction term into a non-interacting 
term coupled to an auxiliary field. This transformation 
enables us to treat Hubbard models with non-local elec¬ 
tron interactions, provided that the long-range interac¬ 
tion gives rise to a transformation matrix that is posi¬ 
tive definite (a non-positive definite transformation ma¬ 
trix corresponds to a diverging auxiliary field). 

In particular, we find that the transformation matrix 
for the case where the Vij is obtained from the QC-PPP 
method is not positive definite. This is a direct con¬ 
sequence of the interpolation scheme mentioned above, 
which renders the off-diagonal matrix elements associ¬ 
ated with the non-local interaction comparable with the 
diagonal elements associated with the local interactions. 
As a result, we could not use quantum Monte Carlo to 
simulate the QC-PPP model. 

In the DQMC, we used inverse temperature /3 = ^ in 
Eq. [5] between 24 and 36, which is sufficient to probe the 
low-energy behavior of the system. In the PQMC, we 
chose the Hartree-Fock state as our trial wave-function, 

| Tt), using 0£ = 40 (see Eq. [6| to project the wave- 
function onto the ground state. We made use of the 
scaling behavior of the antiferromagnetic structure factor 
( Safm ) to estimate the magnetic state of the system. 


Safm J 2 ^{ m iArrijA} + ( miBrrijB) 


(7) 
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FIG. 3. Antiferromagnetic (AFM) order parameter rriAFM — 
V Safm /(2L 2 ) (on top) and scaled antiferromagnetic struc¬ 
ture factor SafmL 2 ^/ 1 ' /N (bottom) in terms of the inverse 
system size. We have used both projective quantum Monte 
Carlo (blue) and determinant quantum Monte Carlo (red) to 
study the phase of graphene subject to: 18% biaxial strain 
within the Thomas-Fermi model (upper points); 18% biax¬ 
ial strain within the constrained RPA model (middle points); 
0% strain within both the Thomas-Fermi and the constrained 
RPA models (lower points). We could not simulate the quan¬ 
tum chemistry - Pariser-Parr-Pople model with auxiliary field 
quantum Monte Carlo since its partially screened Coulomb 
potential gives rise to diverging auxiliary fields. 


where rriic stands for the magnetization of the site lo¬ 
cated in the atom of sub-lattice C = A, B of the unit cell 
iq, while L 2 is the number of unit cells (i.e. N = 2 L 2 
sites). The system’s AFM order parameter is given by 
ttiafm = V Safm /(2T 2 )• We have simulated lattice 
sizes between L = 6 and L = 15. In order to take finite 
size effects into account, we utilize ttiafm = m afmL» 
where we use the critical exponents /3/v ~ 0.9 (obtained 
from the best data collapse in Ref. m ), compatible with 
the Gross-Neveu universality class Ham]. 

Figure [ 3 ] shows the dependence of ttiafm and Safm 
with the system size. For unstrained graphene both 
the cRPA and TF methods show Safm decreasing with 
system size (and ttiafm extrapolating to zero in the 
thermodynamic limit L -A 00 ), indicating that with¬ 
out strain suspended graphene is metallic (in agreement 
with experimental observations). However, most inter¬ 
estingly, with 18% biaxial strain, the TF model shows 
Safm increasing with the system size (with ttiafm ex¬ 
trapolating to a non-zero value when L -A 00 ), in¬ 
dicating an anti-ferromagnetic Mott insulator in the 
thermodynamic limit. This corresponds to an interac- 
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tion driven gap of A = (0.55 =b 0.05) eV, compara¬ 
ble to estimates in Ref. m obtained by hybrid den¬ 
sity functional calculations (Hartree-Fock exchange hy¬ 
bridized with generalized gradient approximation for the 
exchange-correlation) that do not accurately treat the 
effects of electron correlations. Moreover, within the 
Thomas-Fermi method, our QMC calculations find a 
critical strain of rj c « 0.15. Notice that in this case, 
U — U — V — 3At < U c demonstrating that the sug¬ 
gestion by Ref. m for mapping the long range Hub¬ 
bard model for graphene into an effective onsite Hubbard 
model is quantitatively inaccurate. 

Although the TF method has no adjustable parame¬ 
ters, it assumes that the Coulomb interaction between 
p z electrons on the same atom and between neighboring 
atoms is screened in the same way [26]. This assumption 
slightly overestimates the ratio U/V giving a smaller crit¬ 
ical strain r] c for the SM-AFM transition. On the other 
hand, the canonical cRPA method ignores bandwidth 
and low-energy spectral weight reduction originated from 
integrating out the high-energy bands [32]. This gives 
rise to artificially weak partially screened Coulomb in¬ 
teractions, resulting in an overestimation of the critical 
strain r\ c . Due to finite sizes, the PPP model overesti¬ 
mates the value of U and V, and the Ohno interpolation 
underestimates their difference. However, looking at the 
three models together, we therefore conclude that the 
profile of the Coulomb potential for realistic graphene 
lies somewhere in between the TF and cRPA estimates. 
The TF model gives a maximum Mott gap of more than 
an order of magnitude larger than room temperature, 
and this value sets the upper bound for experiments. 

In summary, using the best available models in the lit¬ 
erature to estimate the effective Coulomb interaction be¬ 
tween pz electrons in graphene, we have employed quan¬ 
tum Monte Carlo simulations to explore graphene’s phase 
diagram in response to parameters that can be changed 
experimentally. We have found, surprisingly, that manip¬ 
ulating the short-range part of the effective Coulomb po¬ 
tential (i.e. the non-universal and material specific com¬ 
ponent) is the crucial factor in determining the phase of 
graphene. Most importantly, we show that application 
of experimentally realistic amounts of isotropic strain is 
a promising route to cross the SM-AFM quantum phase 
transition and to observe a strongly correlated state in 
this otherwise weakly interacting material. 
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